First-Principles Calculations at Constant Polarization 



in 
o 
o 

(N 

> 
O 

(N 



C/3 



X3 

o 



> 



o 

s 

I 

o 
o 



X 



Oswaldo Dieguez and David Vanderbilt 
Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854-8019, USA 

(Dated: February 2, 2008) 

We develop an exact formalism for performing first-principles calculations for insulators at fixed 
electric polarization. As shown by Sai, Rabe, and Vanderbilt (SRV) [N. Sai, K. M. Rabe, and 
D. Vanderbilt, Phys. Rev. B 66, f04f08 (2002)], who designed an approximate method to tackle 
the same problem, such an approach allows one to map out the energy landscape as a function 
of polarization, providing a powerful tool for the theoretical investigation of polar materials. We 
apply our method to a system in which the ionic contribution to the polarization dominates (a 
broken-inversion-symmetry perovskite), one in which this is not the case (a III-V semiconductor), 
and one in which an additional degree of freedom plays an important role (a ferroelectric phase 
of KNO3). We find that while the SRV method gives rather accurate results in the first case, the 
present approach provides important improvements to the physical description in the latter cases. 

PACS numbers: PACS: 71.15.-m, 77.22.-d, 77.80.-e. 



In 1993 King-Smith and Vanderbilt introduced a 
theory for computing the electric polarization of an infi- 
nite solid, showing for the first time that this property 
was indeed a well-defined quantity for an insulating ma- 
terial. Their theory of the bulk polarization (TBP) is 
based on computing the electronic contribution to the 
polarization as a Berry phase of the valence-band Bloch 
wave functions transported across the Brillouin zone. A 
practical numerical scheme to compute it was given in 
the same paper 0], and it is now widely used in first- 
principles calculations. Later, Souza, Iniguez, and Van- 
derbilt (SIV) used the TBP as the cornerstone for a 
method to calculate the exact ground state of an insu- 
lator in the presence of an electric field by minimizing 
an electric enthalpy functional expressed in terms of oc- 
cupied Bloch-like states on a uniform grid of points in 
reciprocal space. It is therefore possible to compute from 
first-principles many interesting properties of materials 
related to their behavior under electric fields. 

In this Letter we propose a method to do first- 
principles calculations not at constant electric field, but 
at constant electric polarization. In this way, it would 
become possible to map the energy E of an insulator 
as a function of its polarization P. There are several 
reasons why this is useful. First, it allows for a exhaus- 
tive search for competing local minima and for saddle 
points in a system with a complicated energy surface. 
These features would otherwise be hidden in a standard 
optimization procedure, which would only be likely to 
find a single miminum of the energy. Second, know- 
ing E{P), we can calculate various properties related to 
derivatives of E with respect to P, such as the linear 
and non- linear dielectric susceptibilities. While methods 
of density-functional perturbation theory (DFPT) Q can 
also be used to access these susceptibilities, our approach 
is advantageous in that it does not require much special- 
purpose programming and yields non-linear susceptibil- 
ities with very little extra effort beyond that needed to 



obtain the linear ones. And third, Landau-Devonshire 
theories have historically provided an important avenue 
to the understanding of ferroelectric materials, based on 
an expansion of the energy in powers of polarization. The 
ability to compute E{P) may open the door to the first- 
principles derivation of Landau-Devonshire descriptions, 
as opposed to the usual empirical formulations, for a wide 
range of ferroelectric materials. 

The problem of performing simulations at constant 
polarization in the context of first-principles calcula- 
tions has been addressed previously. In particular, Sai, 
Rabe, and Vanderbilt (SRV) developed an approxi- 
mate method to do this and applied it to obtain E(P) for 
several perovskite systems of interest. A similar approach 
was used by Fu and Bellaiche to study the electromechan- 
ical response of solids under finite electric fields |S] . Our 
method relies on many of the ideas of the SRV one, but 
it is instead an exact method because it incorporates the 
SIV theory of finite electric fields . 

Our goal is to find the minimum energy of a sys- 
tem of atoms for which the electric polarization takes 
some target value Pt. The system has both ionic and 
electronic degrees of freedom, and we will use vectors X 
and ^ to represent them. Vector X has 3iV + 6 compo- 
nents given by the cartesian coordinates of the N atoms 



FIG. 1: Sketch of energy as a function of polarization for a 
double-well potential. The value of the electric field needed 
to impose the polarization is the same at the three indicated 
points, as their slopes are the same. 
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in the simulation cell and the 6 components of the strain 
tensor. On the other hand, ^ contains the one-electron 
wave functions of the system. We are therefore facing a 
standard constrained optimization problem that can be 
solved by introducing a Lagrange multiplier A and search- 
ing for the minimum of E + \ ■ (P — Pt), or, discarding 
the constant term, the minimum of -I- A • P. Defining 
A = — ftS, where is the volume of the unit cell of the 
material under consideration, we have 



min {£;(X,*) + A-P(X,*)} 



min {min{i;(X,*) - 05 • P(X,*)}} 

£^P = Pt 

min H{X,S), (1) 



where A P = Pt indicates a minimization over the 
values of A that impose the constraint P = Pt . The func- 
tion -ff (X, £) is the electric enthalpy that is minimized 
in the SIV method to find the electronic ground-state of 
an insulator in the presence of an electric field £ 

Therefore, the problem of finding the minimum energy 
of a system that is constrained to have polarization P is 
equivalent to finding the ground state of a system under 
an electric field £ that imposes polarization P. However, 
the apparently straightforward approach of doing calcu- 
lations at different values of the electric field to find E{P) 
parametrically will not work in general. Fig. ^illustrates 
why this is so for a typical double-well potential charac- 
teristic of ferroelectric perovskites. In this case, several 
values of the polarization correspond to the same value 
of the electric field, since £ — n^-^{dE/dP), and an op- 
timization using the SIV theory will most likely find the 
point that has the lowest energy of the three, i.e., the 
global minimum of H{P). With extra care it might be 
possible to find the secondary local minimum of H{P), 
but not the third solution, which is a maximum (or, in 
3D, a saddle point) of H{P). It follows from this line of 
reasoning that one cannot map the region of E{P) that 
has negative curvature, and it may also prove difficult to 
map nearby regions of positive curvature. 

We now describe how to perform a minimization of the 
enthalpy H over the ionic degrees of freedom X and the 
electric fields £ that produce a desired polarization Pt 
in a way that parallels the presentation in Sec. III. A of 
Ref.S From now on, we will assume that the optimiza- 
tion is done while keeping the cell vectors fixed. (It is 
possible to remove this constraint, but it is necessary to 
take into account some technical subtleties related to the 
way in which the stress tensor is computed in the presence 
of an electric field. The details, together with examples, 
will be presented elsewhere Q.) We begin with a trial 
guess (Xo,5o), and we expand iJ(X,5) and P(X,£) as 



low-order Taylor series 



H = iJo - F5X+ -(5XK5X- OP(55, (2) 



(3) 



where Hq = ff(Xo,5o) and Pq = P(Xo,5o) (atomic 
units are used throughout unless stated otherwise). The 
Fi = —{dH/dXi) are the 3iV components of the forces 
on the atoms at finite electric field the Kij — 
{d^H/dXidXj) are the 3N x 3iV force-constant matrix 
elements, the Zia — fl{dPa/dXi) are the 3A^ x 3 Born 
effective charges, and the Xap = 4:Tr(dPp/d£a) are the 
3x3 elements of the dielectric susceptibility tensor. We 
can then predict that dH/dXi = and P = Pt at an 
X = Xq + (5X and slti £ = £q + 6£ given by 
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where APq = Pq — Pt- This system of linear equa- 
tions can be solved iteratively, refining (SX and 6£ until 
F and APq both vanish. As mentioned before, F can 
be computed exactly in the presence of an electric field 
according to the SIV prescription 0. The guiding ten- 
sors K, Z and x do not need to be computed exactly; 
replacing these by approximate versions only results in a 
slower convergence to the correct solution, without shift- 
ing the solution itself. In particular, we normally find it 
sufficient to compute K, Z and x at zero electric field. 

The scheme described here has been implemented 
to work with the ABINIT "T] density-functional theory 
(DFT) 1^ code. Instead of performing a direct iterative 
solution of Eq. Q), we have used a nested-loop algorithm 
for the sake of robustness. Starting with some guess for 
X and £, we keep the atoms fixed and vary the field in 
the internal loop until the polarization is the target one, 
a problem that is well behaved. Once this is achieved, we 
solve Eq. Q to get new values of X and £, and iterate 
until convergence is achieved. 

As a first example of how our method works, we apply 
it to a soft-mode system studied by Sai, Rabe, and Van- 
derbilt 01 ■ They have shown that breaking the inversion 
symmetry in perovskites by modulating their composi- 
tion in a cyclic sequence of layers produces some interest- 
ing features in their energy landscape, apart from making 
them promising candidates for new materials with useful 
piezoelectric and dielectric properties. They presented 
results for Ba(Ti-(5, Ti, Ti-|-(5)03, where the two species 
that alternate with Ti on the B site are virtual atoms 
that differ from Ti in that their nucleus has an defect or 
excess of charge equal to 5. For different values of 5 they 
find the E{P) curves that we reproduce in Fig. |21 (dashed 
lines). As i5 increases, the double well potential becomes 
asymmetric (a feature not seen in normal perovskites) 
until one of the minima eventually disappears at around 
6 = 0.4. 
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FIG. 2: Energy versus polarization along the c axis for the 
broken-inversion-symmetry perovskite described in the text, 
for different values of &. The solid (dashed) line is a polyno- 
mial fit of the values obtained using our (the SRV) method. 



We have repeated their study using the same DFT 
methodology, but applying our exact method to com- 
pute E{P) instead their approximate one. We used 
the ABINIT j^] code to do our calculations, with the 
Ceperley- Alder form of the local-density approxima- 
tion lCr| (LDA) to obtain the exchange-correlation term 
in DFT, a plane-wave cutoff of 35 Ha to define the ba- 
sis set, a 4 X 4 X 4 Monkhorst-Pack grid for com- 
putations in reciprocal space, and TrouUier-Martins jl^ 
norm-conserving pseudopotentials (13i] to model the ion- 
electron interactions. The tetragonal unit cell vectors 
were kept fixed, with cell parameters a = 7.547 a.u. and 
cja — 3.036. Here, our results (solid lines in Fig. I^J are 
very similar to the SRV ones; the curves only differ no- 
ticeably for large values of the polarization. This is not 
surprising, since the dielectric behavior is dominated by 
the ionic response for soft-mode systems like this one. 

On the other hand, the electronic response is found 
to be profoundly important to the dielectric behavior in 
the case of a HTV semiconductor like AlAs, as can be 
seen m Fig. 13 Here we employed ABINIT using the 
LDA 1^, a plane- wave cutoff of 9 Ha, a 6 x 6 x 6 recip- 
rocal space grid , TrouUier-Martins |0| pseudopoten- 
tials [lj| , and a theoretically optimized lattice constant of 
a = 10.62 a.u. As can be seen in Fig.O there is a drastic 
difference between the E versus P curves computed using 
the SRV method (ionic response only) and the one found 
with our new method (ionic and electronic response). In 
order to quantify this difference, the dielectric constant in 
both cases was computed from the curvature at the min- 
imum of each function as e = 1 -I- Att / {cf E / dP^) . When 
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FIG. 3: Energy versus magnitude of polarization in the line 
joining two nearest Al and As neighbors. The solid (dashed) 
line is a polynomial fit of the values obtained using our (the 
SRV) method. 



computed with the SRV method, e = 3.0, while when we 
use our new method, e = 10.3. The experimental result 
is e = 10.1 0. 

Our last example involves potassium nitrate, which has 
an interesting phase diagram that includes a reentrant 
ferroelectric phase (phase HI, RSm), and that has been 
proposed as a promising material to be used in random- 
access memory devices 15]. The unit cell of phase HI 
KNO3 is a five-atom rhombohedral cell that results in 
alternating planes of K atoms and NO3 groups arranged 
in such a way that the K atom is not equidistant from the 
NO3 groups above and below it, but instead is slightly 
displaced in the vertical direction, giving rise to a polar 
structure. Figs. 01a)-(b) show the conventional 15-atom 
hexagonal cell that is most convenient for visualization. 
The calculations are performed using ABINIT |7(|, the 
LDA , a plane- wave cutoff of 30 Ha, a reciprocal space 
grid with 6 inequivalent points, and TrouUier-Martins 
[12.] pseudopotentials . The theoretical structural op- 
timization of bulk ferroelectric KNO3 gives hexagonal 
lattice parameters of a = 9.68 a.u. and c = 14.86 a.u., 
to be compared with the ones found experimentally at 
91 °C, a = 10.37 a.u. and c = 17.30 a.u. As for 

the internal degrees of freedom, the distance between a 
K atom and the N atom above it is found to be z = 0.57c 
(experimentally, z = 0.59c while the N-O distance 

is d = 2.35 a.u. (experimentally, d = 2.42 a.u. 01 )• The 
computation of the polarization gives 0.16 C/m^, higher 
than the experimental values of 0.08-0.11 C/m^ given in 
Ref. 18. (This disagreement is related to the reduction in 
volume found theoretically; when we performed an opti- 
mization of the atomic positions using the experimental 
lattice constants, we found a polarization value of 0.075 
C/m2.) 

Fig-EJc) shows the energy of the ferroelectric phase of 
KNO3 as a function of the magnitude of the polarization 
along the c axis. As in the case of perovskites, we can see 
a double well potential, but now the change of polariza- 
tion given by the movement of the NO3 atoms parallel to 
the c axis is accompanied by a rotation of these groups in 
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FIG. 4: (a) Side view of the ferroelectric phase of KNO3 
(phase III), with the c hexagonal lattice axis running verti- 
cally, (b) Top view of the same structure, (c) Energy versus 
polarization P — (0,0, P). (d) Angle of rotation of the NO3 
group as a function of polarization P. 



their plane, as shown in Fig. ^d). When doing SRV cal- 
culations, both the derivative of the energy and the angle 
of rotation are continuous functions of the polarization, 
and the NO3 group rotates from 0° to 60° as the polariza- 
tion goes from negative to positive, with the paraelectric 
configuration that corresponds to the maximum of E[P) 
being the highly symmetrical one in which z = 0.5 and 
the rotation angle is 30°. However, when the new exact 
method is used, this behavior changes noticeably: the 
E{P) curve has a discontinuous derivative at P = 0, and 
the angle-rotation curve is no longer continuous. 

Although it may seem puzzling that the response of 
the system changes so significantly just by including 
the electronic response, this behavior can be understood 
on the basis of a simple model in which the energy is 
a function not only of polarization P, but also of the 
rotation angle 6. To understand the qualitative fea- 
tures, it is sufficient to consider a low-order expansion 
E{P, 9) = cos W + a cos 129 + (3 P cos 39 + P^. We as- 
sume that P > 2-\/2, in which case E{P) has two minima. 
Then, depending on the value of a, the behavior can be 
continuous {a < 1/4) or discontinuous {a > 1/4) in P. 
Here, it appears that the system was already close to 
this critical value of a, and inclusion of electronic effects 
happened to shift the system from the continuous to the 
discontinuous regime. A more detailed discussion of this 
material, and its modeling along the lines sketched above, 
will be given elsewhere. In any case, the ability of our 
approach to describe the complexity of the structural be- 
havior of KNO3 under polarization reversal provides an 
excellent example of the power of the method. 

To summarize, we have presented a method for find- 
ing the most stable structural configuration of an insulat- 



ing crystal when its electric polarization is constrained to 
take on a given value. Our method builds upon an earlier 
approach [J| that makes the approximation of including 
only the lattice response to applied electric fields. By us- 
ing the recently-developed theory of finite electric fields 
to include also the electronic response of the system, 
we have developed a new approach that is, instead, ex- 
act. We have illustrated the method by obtaining E{P) 
curves for three rather different kinds of insulating ma- 
terials, and have illustrated how the method is capable 
of describing the complexity of the nonlinear structural 
and dielectric response. 
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